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An approach to compute the anharmonic peaks of the phonon dispersion curves through the 
ah initio calculated Hellmann-Feynman forces from a series of supercells with realistic atomic dis¬ 
placements of all atoms, which correspond to a given temperature, is reported. Obtained phonon 
dispersion bands are able to represent the positions and shapes of the anharmonic peaks. As exam¬ 
ple, the approach to cubic PbTe and perovskite MgSiOs crystals is applied. 
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The conventional treatment of anharmonicity in crys¬ 
tals relies on the expansion of the potential energy with 
respect to atomic displacements up to orders higher then 
quadratic terms. Typically, third and fourth orders ex¬ 
pansion terms are used [1], in spite of existence of pow¬ 
erful perturbation methods being able to sum part of 
expansion diagrams mil]. Unfortunately, these methods 
are computational very expensive. 

The anharmonicity is required for derivation of the 
thermal conductivity of materials, specially thermo¬ 
electrics, which efficiency to transfer heat to electricity 
increases with lowering the heat transfer. Moreover, the 
high pressures thermal conductivity is important to un¬ 
derstand heat transfer in the Earth interior. 

In this Letter a new approach to handle crystals anhar¬ 
monicity is proposed. It relies in probing the potential 
energy landscape of the crystal with atomic vibration. 
Then, the vibrations are transformed and classified in 
the wavevector reciprocal space, similar to harmonic case. 
The expansion of the potential energy is not needed. 

The existing density functional theory (DFT) codes 
provide sufficiently accurate Hellmann-Feynman (HF) 
forces to calculate the phonon frequencies and properly 
describe related phenomena. The DFT approach requires 
to approximate the crystal as a supercell, which is a par¬ 
allelepiped construction made of primitive unit cells. By 
definition periodic boundary conditions are imposed on 
the supercell. 

To compute harmonic phonons, the ab initio force con¬ 
stant approach formulated by Parlinski, Li, Kawazoe, in 
1997 [5], can be used. This method has been already ap¬ 
plied for several hundreds of crystals and the results show 
good agreement with the measured phonon data. The 
computing method relies in using the relation between 
harmonic forces FH(n,/i), induced by atoms displaced 
by u(m, u) from the equilibrium position 

FH(n,/u) =• u(m,i/) (1) 

where n, m label the primitive unit cells in 
the supercell, and fi, v number atoms within 
these cells. The BH(n,/i, m, z/) are harmonic 


force constants between marked atoms and they 
determine the conventional dynamical matrix 
= -^^^^^BH(0,Ai,m,Oea;p{-27rk • 

{R-(o,m) - R.(m,0]} El- The eigenvalue equation 
(k, j)£(k, j) = D(k)f(k,ji) of this Hermitian matrix 
provides phonon frequencies and eigenvectors 

f(o))(k, j). In Ref. [5] was proposed to decouple each 
force constant matrix into a product 

BH(n, /i; m, u) = A(n, fi; m, n) • PH(n, /i; m, u) (2) 

where (3 x 3) BH(n,/i, m, z^) matrix is rearranged to 
(9 X 1) matrix, A(n, //; m, z/) is a (9 x p) matrix, and 
PH(n, /i; m, z/) is a (p x 1) matrix. The A(n, p; m, z/) 
is entirely determined by crystal symmetry, and does 
not involve potential parameters. The parameter matrix 
PH(n, p; m, u) depends on potential details only, and in¬ 
dex p denotes number of independent parameters of the 
force constant in question. Rearranging the three compo¬ 
nents of displacement vector u(m, u) into (3 x 9) matrix 
U(m, u), Eq. 0 becomes 

FH(n,M)= F Cu(n,/i,m,iy) ■ PH(n,/i;m,iy) (3) 

m,u,j 

where the (9 x p) symmetry adapted displacement matrix 
Cu(n, /i, m, u) has been determined as 

Cu(n,/i;m, z/) = -U(m, z/) • A(n, p; m, z/). (4) 

The Cu(n, p; m, z^) matrix is known from used atomic 
displacements and crystal symmetry. For small atomic 
displacements of order of 0.03A, the harmonic forces 
FH(n,/i) can be obtained from ab initio calculated HF 
forces. These data permit to find all unknown indepen¬ 
dent parameters PH(n, p; m, z^), and from Eq. ([^, all 
force constants BH(n, p; m, u) within the supercell. 

Single atomic displacements in the primitive unit cell 
(0, u) is sufficient to create single HF force field for that 
atom. A force field is a set of 3n HF forces obtained for 
the run in ab initio code with a single displaced atom. 
Here, n denotes the number of atoms in the supercell. 
Usually a single HF force field is not sufficient. A mini¬ 
mum number s of HF force fields is equal to the number 
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of non-equivalent atoms of primitive unit cell, supple¬ 
mented by a number of non-equivalent directions of the 
displacements. Finally we have 3ns equation of kind @ 
to be solved simultaneously. Collecting these 3ns equa¬ 
tions to a global form, the system of equations for the 
harmonic force constant parameters can be written as 

= Cu • Vn (5) 

where Cu and Vn are (3ns x 1), (3ns x p') and 
(p' X 1) dimensional matrices, respectively, where p' is 
now the total number of independent parameters needed 
for all force constants within a super cell. In this sys¬ 
tem of equations the number of HF forces is greater then 
the number of potential parameters, 3ns > p', therefore, 
this is an overdetermined system. To solve it, the sin¬ 
gular value decomposition method (SVD) [6] to matrix 
Cu is applied [5]. Then the potential parameters can 
be found from HF force Vn = Cu~^ • J^n- This SVD 
method provides a solution, which is the best approxi¬ 
mation in the least square sense. The above procedure 
delivers very effective method to find harmonic phonons, 
their frequencies, polarization vectors, phonon dispersion 
curves, phonon density of states, and many other phonon 
dependent quantities. 

The global form of system of equations, Eq. <§ , per¬ 
mits to add other conditions to be fulfilled. The simple 
example are the translational-rotational invariants [3], 
which can be reformulated to a form of matrix Ad of 
(18n X p') dimensions using Eq. (H- Then, Eq.0 be- 
comes 


where (3 is adjusting the strength to satisfy the 
translational-rotational conditions. Eq. is solved by 
SVD method. 

To compute anharmonic effects one may use a similar 
approach. Any DTE calculations of HE forces contain 
an information on the anharmonicity. This means that 
there is an access to anharmonic landscape of the crystal 
potential energy. Generally, a HE force F(n,/i) acting 
on an atom (n, p) can be treated as originating either 
from harmonic FH(n,/r), or from anharmonic FA(n,/i) 
contributions. In the pure harmonic regime the matrices 
BH(n,/i; m, z/), Eq.Q, do not depend on the amplitude 
of atomic displacements, hence they are temperature in¬ 
dependent. It follows from linear dependence between 
harmonic forces and atom displacements. 

The anharmonic potential has a many-body character, 
involving at a given temperature a lot of atoms simulta¬ 
neously displaced. Additionally, all atoms of the crystal 
move in time, varying the displacements around the equi¬ 
librium position, and at a given time moment all atoms 
of the supercell create pattern of displacement de¬ 
noted as (i). One may treat also the displacement pattern 


(i) as coming at the same time from different locations 
of the crystal. Each displacement pattern (i) allows to 
calculate one HE anharmonic force field F^\ Knowing 
the displacement patterns and force fields F^\ we 
supplement the Eq. by anharmonic contributions, for¬ 
mulating in this way a global anharmonie form of system 
of equations 



\ ( Cu \ 
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where i = 1, 2,... V runs over displacement patterns. 

The SVD solution of Eq. 0, corresponding to single 
displacement pattern, gives a single set of dispersion 
curves of index (i). Since in space, or in time the dis¬ 
placement pattern changes, for a realistic modelling 
of anharmonic effects one should select a sequence of N 
such patterns (i = 1,2, ...V), and compute dispersion 
curves for all these N patterns. Erom set to set, prob¬ 
ing different anharmonic environments, the phonon dis¬ 
persion curves are slightly different due to variation of 
potential energy reached now by displaced atoms. In¬ 
deed, in the absence of anharmonicity, for example with 
negligible amplitudes of atomic displacements, all sets (i) 
of phonon dispersion curves will be identical to harmonic 
case within the accuracy of a computational noise and ap¬ 
proximate SVD solution. In anharmonic case all phonon 
dispersion curves, i = 1, 2,... V form bands, which re¬ 
flect the anharmonic character of the potential energy 
landscape, caused by the atomic displacements around 
the global potential minima. 

To stay consistent with the supercell concept, the 
displacement patterns must preserve the periodic 
boundary conditions. Thus, a displacement pattern must 
be a superposition of the allowed phonon displacement 
waves. Hence, each phonon displacement wave must be 
commensurate with the supercell size and shape. Its am¬ 
plitude can be estimated from the harmonic theory as 


U(n,/x) 


e(k,j) 

yMp 


exp[27ri('k • R(n, p) — 0(k, j)] 


(8) 


where the mean square amplitude of the wave is deter¬ 
mined by 


< Q^(k,i) >= 


h 

2w(k,ji') 


coth 


/ fko{k,j) \ 

V ^keT ) 


(9) 


and are the harmonic phonon frequencies at 

wavevector k, and phonon branch j. The phase factor 
0(k,jf) could be taken at random to mimic different dis¬ 
placement patterns. The temperature displacement vari¬ 
ation, Eq. 0 depends on the thermal occupation factor, 
Eq. 0- The HE forces calculated for a displacement pat¬ 
tern include contributions arising from many atoms dis¬ 
placed simultaneously. The approximation relies on us¬ 
ing only displacement waves of a few hundreds discreet 
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wavevectors. We remark that the displacement patterns 
may also be obtained in another way, for example, as 
snapshots of atomic motion traced during molecular dy¬ 
namic simulation. But even in this case the confinements 
imposed by selecting only commensurate displacement 
waves holds as well. 

The crystal symmetry is determined by the crystallo¬ 
graphic space group. In harmonic regime all fluctuations 
- atomic displacements - are described and classified by 
the phonon modes. To each phonon mode an irreducible 
representation of the crystal space group is assigned. 
The above classification applies to phonon modes derived 
from force constants being a solution of Eq.J^ , and each 
set (i) of global system of equations, Eq.([^. Then, it 
follows that the current anharmonic phonon bands are 
characterized by the same irreducible representations as 
the irreducible representations of corresponding phonon 
dispersion relations of the harmonic case. However, the 
anharmonic phonon bands have much larger chances to 
overlap. Therefore, a systematic method to select out a 
single anharmonic peak from phonon bands is required. 

Now the projection method to select out a particu¬ 
lar anharmonic mode is proposed. The conventional 
assignment of irreducible representations to harmonic 
phonons is considered as a reference to classify the an¬ 
harmonic peaks. Then, diagonalization of the dynamical 
matrix, delivers orthonormalized eigenvectors f J), 

and of harmonic phonons curves and anhar¬ 

monic phonon bands relevant for all displacement pat¬ 
terns (i), where J and j label phonon modes for the same 
wave vector k, respectively. Each eigenvector involved in 
the anharmonic peak can be expanded over complete set 
J = 1,2,... Jmax of harmonic eigenvectors 

Jmax 

^ a«(k,j, J)f(°)(k, J). (10) 

J=1 


Applying the orthonormality relation f*^*^(k,jf) • 

f^*^(k,ji) = 1, the expansion coefficients of Eq.(lO) can 
be found as (a^*^(k, j, J) = f *^*^(k, j) ^^^(k, J). At fixed 
k wavevector the anharmonic J mode of phonon density 
of states^ denoted by bj{uj, k), as a function of frequency 
cj, can be found from the histogram 


-J N jmax 

I J) p (5Au;(w-w(*))(k,i)) 

i=l j=l 

( 11 ) 

where Z = N • Jmax 'jmax ' The histogram bin Acj is 
defined by the function Sauj(^) = 1, T — ^ < x < or 
0 otherwise. The summation i = 1, 2,... A runs over all 
displacement patterns (i). The coefficients (a^*^(k,j, J) 
select out from all (i) bands only those phonon inten¬ 
sities which resemble the vibrations determined by har¬ 
monic eigenvectors J). Each anharmonic mode 

of phonon density of states bj{uj,'k) determines a single 





FIG. 1. (Color online) The cubic PbTe phonon dispersion re¬ 
lations along r(0, 0, 0) — A(|, 0, 0) wavevector path, (a) Har¬ 
monic curves at T = OK (temperatureless regime), and pres¬ 
sure P = OGPa. (b) Anharmonic phonon bands originated 
from 100 different sets of displacement patterns at T = 300A 
and P = 0.828GPa. (c) Anharmonic mode density of states 
of all phonons, Eq( 11), for a fixed wavevector k = A (vertical 
z-z line at (b) plot), and as a function of phonon freq uency u) 
in THz. The total spectrum was separated, Eq. 1^, to TA, 
LA, TO and LA anharmonic peaks. 


anharmonic mode of symmetry J. The referred method 
is able to select out all single anharmonic peaks, even if 
they overlap. 

The above approach is illustrated with two examples: 
PbTe and MgSiOs crystals. Our intention, however, is 
not to provide a deep physical analysis of the samples, 
but to show what kind of results can be obtained. All nu¬ 
merical calculations were performed within DET method 
using VASP code [7], and applying GGA-PAW approach 
issued with this code. Phonons were calculated with 
PHONON code [5]. Displacements of 0.03 A were used 
to generate the lists of harmonic HE forces. 

The anharmonicity of cubic PbTe semiconductor has 
been extensively studied in Ref. [8]. Eor PbTe, the space 
group Fm3m, 2x2x2 supercell with 64 atoms, and 
4x4x4 /c-point mesh were used. The optimized lattice 
constant for pressure P = OGPa was a = 6.557A. The 
effective charges of Z* = ±5.8 and electronic dielectric 
constant e = 25.26 were taken from Ref. [9]. The har¬ 
monic phonon dispersion relations are shown on Eig{^ 
along r(0,0,0) — A(|,0,0) wavevector path. The com¬ 
puted eigenvectors specify uniquely transverse (T), longi¬ 
tudinal (L), acoustic (A), and optical (O) polarizations. 
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indicated as TA, LA, TO, LO modes on the plot. 

For temperature T = SOOLC, 100 displacement patterns 
were generated using in each case 189 phonon waves with 
wavelengths commensurate with 2x2x2 supercell. For 
each pattern the phonon phase was taken at random. The 
average atomic displacements over all patterns were 0.16 
and 0.12 A, and from analytical Eq.M) with harmonic 
density of states were 0.15 and 0.12 A, for Pb and Te, 
respectively. Displaced atoms increase the potential en¬ 
ergy of the pattern with respect to T = OLC state by AE. 
Average value of AE, calculated from the harmonic form 
of the potential energy derived using the same force con¬ 
stant as involved in phonon bands, and on other hand 
from the energy excess provided by VASP, were pretty 
close, and read 75.37 and 74.13 meV/prim.unit cell, re¬ 
spectively. Moreover, the displaced patterns, according 
to VASP, created a pressure of 0.828 GPa in the super¬ 
cell with the same lattice constant as in T = OLC, namely 
a = 6 . 557 A . It is natural since the calculation at T = OLC 
and T = 300/c were done at the same supercell volume. 

Fig|^ shows plots of phonon bands calculated for 100 
displacement patterns. Each band consists of 100 phonon 
curves (200 curves for doubly degenerate modes) of the 
same symmetry. There are regions were bands over¬ 
lap. Having 600 modes for each wavevector k one may 
construct a histogram of the anharmonic mode phonon 
density of states. Such a histogram, plotted for a sin¬ 
gle wavevector k = V(^,0,0) (line z-z), is presented on 
Eigj^. It shows total mode phonon density of states and 
its separation into four specific anharmonic modes TA, 
LA, TO, and LO. The TA and LA modes look rather 
sharp, while TO and LO are much wider. Since the shape 
of each anharmonic peak is known, there is no problem 
to establish peak position and width. Without any addi¬ 
tional ah initio calculations, but for the same T and P, 
the anharmonic peaks for any wavevector k can be found. 
One notices that on Eigj^ imaginary (negative) phonon 
branches appear, however, they are very rare, and their 
contributions to anharmonic mode of phonon density of 
states, FigjlJ;, is negligible. 

The anharmonicity of MgSiOs has already been car¬ 
ried on with respect to thermodynamical properties and 
thermal conductivity m of the Earth lower mantle, 
where high T and P persist. The perovskite MgSiOs be¬ 
longs to orthorhombic space group Pmnh^ No=62. Eor 
all runs the 1 x ^/2 x \/2 super cell, and 2 x 2 x 2 k- 
point mesh were selected. Eor effective charges the for¬ 
mal charges were chosen. The anharmonic calculations 
started from optimization of the supercell at pressure of 
P = bl.SGPa and T = OK. Hence, the lattice constants 
a = 6.481A, b = 4.689A, c = 4.462A were found. The 
harmonic phonon dispersion relations were plotted along 
X(|, 0,0)—r(0, 0, 0)—V(0, 0) wavevector path, and are 

shown on Eigj^. Since the primitive unit cell contains 
20 atoms, therefore there appears 60 phonon dispersion 
curves. 
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FIG. 2. (Color online) The orthorhombic perovskite MgSiOs 
phonon dispersion curves along A(|,0,0) — r(0,0,0) — 
y(0, |,0) wavevector path, (a) Harmonic dispersion rela¬ 
tions at T = OK (temperatureless regime), and pressure 
P = 57.3GPa. (b) Anharmonic phonon bands originated 

from 100 different sets of displacement patterns at T = 2300 A 
and P = 70.5GPa. (c) Anharmonic mode density of states 
of all phonons, Eq(pT|, for a fixed wavevector close to k = 
r(0,0,0)) (vertical z-z line at (b) plot), as a function of 
phonon frequency uj in THz. (d) Plot of all anharmonic peaks 
of Big symmetry at F selected out using Eq.(ll), from com¬ 
plete anharmonic mode density of states (line z-z at (b))plot) 

at k = r(0, 0,0). 


For the anharmonic analysis temperature of T = 
2300A was chosen, and 100 displacement patterns was 
generated. Each pattern used 117 phonon waves with 
wavelength commensurate to the supercell size. The av¬ 
erage atomic displacements over all patterns were slightly 
anisotropic, but the average values of them were: 0.125, 
0.105, 0.128, and 0.089 A, and from analytical Eq.(|^) 
with harmonic density of states were 0.110, 0.105, 0.127 
and 0.087 A, for 01, 02, Mg and Si atoms , respectively. 
The potential energy increased due to presence of the 
atomic displacements. The energy excess AE calculated 
directly from force constants and from VASP were 6.049 
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and 6.727 eV/primitive unit cell, respectively. The dis¬ 
placement patters created the pressure P = 70.5GPa for 
the same volume as used in T = OK. 


Next, the anharmonic properties for T = 2300P and 
P = 70.5GPa were calculated. Figj^ shows plots of 
phonon bands calculated along similar wavevector path 
as for harmonic phonon curves. Each band consists of 
100 phonon curves of the same symmetry. There is so 
many bands that they overlap, and the information be¬ 
comes obscure. Moreover, there are observed phonon 
curves which exists in the imaginary frequency region, 
although their intensities were negligible. Having 6000 
modes for each wavevector k one may construct a his¬ 
togram of the total anharmonic mode phonon density of 
states. Such a histogram, plotted for a single wavevec¬ 
tor close to k = r(0,0,0) (line z-z at (b) plot), is pre¬ 
sented on FigfT]^. Applying the projection method, Eqs 
Pini ) to the total anharmonic mode phonon density 
of states, the separation of all 60 anharmonic peaks 
could be done. Eor example, choosing from 60 T-modes 


lAg + ^Big + 5P2^ + P P lOPi^i + IOP 2 U + 
one kind, for example Big and projecting the anharmonic 
eigenvectors onto harmonic eigenvectors, the Big anhar¬ 
monic peaks were found and plotted on Eigj^. All peaks 
are well defined, and could be used to determine their 
positions and widths, and hence the phonon lifetimes. 
Such a procedure of selection anharmonic modes can be 
applied to any wavevector. 


In conclusions we turn attention on several as¬ 
pects: (i) the phonon lifetime r(k, J) is usually found 
from the width of the anharmonic peak r“^(k, J) = 
width^^{hj{uj^\P)^ where operator widths denotes a pro¬ 
cedure to elucidate the width from fitted Gaussian, 
Lorentzian, or similar function. Some of the peaks may 
contain low intensity tails, which influence the phonon 
lifetime, (ii) the ab initio code typically is handling 
phonon-phonon, electron-phonon, and magnetic-phonon 
interactions, and these effects are included in the present 
method, (iii) the anharmonicity is measured by spec¬ 
troscopic methods. The description of the cross sections 


are as a rule supplemented by certain form factor, and 
corresponding theoretical expression can be derived for 
phonon bands. It means that the expected spectroscopic 
spectra can be determined from phonon bands, (iv) for 
the Raman scattering, or infrared absorption the current 
method offer to include changes of Raman tensor itself, 
or effective charges itself due to atomic displacements in 
the patterns, respectively. This effects might influence 
the Raman or infrared spectra, (v) anharmonic displace¬ 
ment patterns offer full parallelization. The CPU time on 
the cluster for MgSiOs example was about 10 hours only. 
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